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Thermally-induced Phases in an Ising Kondo Lattice Model on a Triangular Lattice: 

Partial Disorder and Kosterlitz-Thouless State 
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Magnetic and electronic properties of a Kondo lattice model with Ising localized spins are studied 
on an isotropic triangular lattice. By using Monte Carlo simulation, we present that the model shows 
a rich phase diagram with four dominant states: two-sublattice stripe, three-sublattice ferrimganetic, 
partially disordered, and Kosterlitz-Thouless like quasi-long-range ordered states. Among them, 
the partially disordered state and Kosterlitz-Thouless like state are intermediate phases induced by 
thermal fluctuations in the phase competing regime; they are present only at finite temperatures 
and eventually taken over by another phases as the temperature is further lowered. Although 
the Kosterlitz-Thouless like state was found also in triangular Ising antiferromagnets with further- 
neighbor interactions, the partially disordered state has not been reported in the localized spin 
only models in two dimensions. Interestingly, the partially disordered phase is also peculiar in 
the charge degree of freedom of itinerant electrons; it is insulating and accompanied by charge 
disproportionation. From a combined analysis of a mean-field calculation of the band structure and 
Monte Carlo simulation, we conclude that the partial disorder in the present model is stabilized by 
the Slater mechanism. 

PACS numbers: 75.30.Kz,75.10.-b,75.40.Mg 



I. INTRODUCTION 

The antiferromagnetic (AF) Ising model on a trian- 
gular lattice is one of the most fundamental models for 
geometrically frustrated systems. When the interaction 
is restricted to the nearest-neighbor (NN) pairs, frustra- 
tion in each triangle prevents the system from forming a 
long-range order (LRO) down to zero temperature, and 
the ground state has extensive degeneracy and associ- 
ated residual entropy^ - — . The degenerate ground state 
is extremely sensitive to perturbations. For instance, an 
infinitesimal second-neighbor interaction lifts the degen- 
eracy and induces a LRO in the ground state; a two- 
sublattice stripe order [Fig.QJa)] is selected as the ground 
state when the additional interaction is AF, while a thrce- 
sublatticc fcrrimagnctic (FR) order [Fig. QJb)] is selected 
for the ferromagnetic (FM) interaction. 

In such a degenerate situation, thermal fluctuations 
also play an interesting role. In general, there is a pos- 
sibility that a high-entropic state is selected out of the 
ground state manifold by raising temperature — this is 
called the order by disorder—. For the AF Ising model, a 
candidate for such an emergent state is a partially disor- 
dered (PD) state. The PD state is peculiar coexistence of 
magnetically ordered moments and thermally-fluctuating 
paramagnetic moments. Such possibility was first dis- 
cussed by the mean-field study in the presence of second- 
neighbor FM interaction^; the mean-field study predicted 
that a three-sublattice PD phase with an AF ordering on 
the honeycomb subnetwork and paramagnetic moments 
at the remaining sites [Fig. [TJc)] was induced at finite 
temperature from the degenerate manifold in the limit 
of vanishing second-neighbor interaction. Although such 
PD state was experimentally observed in several Co com- 
pounds^ and theoretically shown to present in a stacked 




FIG. 1. (Color online). Schematic pictures of (a) stripe or- 
der, (b) ferrimagnetic (FR) order, and (c) partial disorder 
(PD) on a triangular lattice. The arrows show magnetically 
ordered sites and the open circles are thermally fluctuating 
paramagnetic sites. 



triangular lattice model^, Monte Carlo (MC) simulations 
in two-dimensional triangular lattice models have indi- 
cated that PD is fragile and remains at most as a quasi- 
LRO; namely, in most cases, the PD state is taken over 
by another peculiar intermediate state, the Kosterlitz- 
Thouless (KT) stateSrii. 

On the other hand, recently, the authors have studied 
Ising-spin Kondo lattice models on a triangular lattice^ 
and kagome lattice^ by MC simulation, and showed 
the presence of PD state in the purely two-dimensional 
models. In these models, the interplay between lo- 
calized moments and itinerant electrons plays a cru- 



cial role in the following points. First, the kinetic mo- 
tion of electrons induces effective interactions known 
as the Rudcrman-Kittcl-Kasuya-Yosida (RKKY) mech- 
anismiZr— . The long-ranged and oscillating nature of 
the interactions drives keen competition between differ- 
ent magnetic states. Furthermore, the change of mag- 
netic states affects the electronic state in a self-consistent 
manner through the spin-charge coupling; the system can 
gain the energy by forming some particular electronic 
state associated with magnetic ordering. In the previous 
study, the authors suggested that the PD state is stabi- 
lized by the non-perturbative role of itinerant electrons^. 

In this contribution, we present our comprehensive nu- 
merical results on the magnetic and electronic proper- 
ties of the Ising-spin Kondo lattice model on a triangu- 
lar lattice. To further clarify the stabilization mecha- 
nism of PD, we analyze the evolution of band structure 
under the PD type magnetic texture on the basis of a 
simple mean-field argument. The analysis suggests that 
the spin-charge coupling can stabilize the PD state by 
the Slater mechanism. Bearing this mean-field picture 
in mind, we present and discuss the results of MC sim- 
ulation in details. We distinguish the two intermediate- 
temperature states, PD and KT-like states, from the two- 
sublattice stripe and three-sublattice FR LRO states, and 
identify the range of the phases by varying the electron 
filling and the strength of spin-charge coupling. Analyz- 
ing the phase diagram and electronic states in compar- 
ison with the mean-field picture, we conclude that the 
two-dimensional PD state is stabilized through the Slater 
mechanism. 

The organization of this paper is as follows. In Sec. [Til 
we introduce the model and method. The definitions 
of physical quantities we calculated are also given. In 
Sec. IIII1 we present the mean-field analyses on the band 
structure in the PD state. MC results are presented for 
magnetic properties in Sec. lIVI and for electronic proper- 
ties in Sec.|Vl Section WT\ is devoted to summary. 



The first term represents hopping of itinerant electrons, 
where Ci a (c ia ) is the annihilation (creation) operator of 
an itinerant electron with spin a =j" , 4- at ith site, and t is 
the transfer integral. The sum (i,j) is taken over nearest- 
neighbor (NN) sites on the triangular lattice. The sec- 
ond term is the onsite interaction between localized spins 
and itinerant electrons, where o\ 



4f c *T 



C \i°i 



resents the z-component of itinerant electron spin, and 
Si = ±1 denotes the localized Ising spin at ith site; J is 
the coupling constant (the sign of J does not matter in 
the present model). Hereafter, we take t — 1 as the unit 
of energy, the lattice constant a = 1, and the Boltzmann 
constant fcs = L 



B. 



Monte Carlo simulation 



To investigate thermodynamic properties of the model 
([T]), we adopted a MC simulation which is widely used 
for similar models^. The model belongs to the class of 
models in which fermions are coupled to classical fields. 
For this class of models, the partition function is given 
by 



Z = Tr f Tr c exp\p(H-nN e )], 



(2) 



where j3 = \/T is the inverse temperature, /x is the chem- 
ical potential, and N e is the total number operator for 
fermions. Here, Tr/ is the trace over classical degree of 
freedom (in the current case, Ising spin configurations), 
and Tr c is the trace over itinerant fermions. In the MC 
simulation, Try is calculated by using the Markov-chain 
MC sampling. MC updates are done by the usual single- 
spin flip on the basis of the standard METROPOLIS 
algorithm. The MC weight is calculated by taking the 
fermion trace Tr c for each configuration of classical vari- 
ables in the following form, 



P({SJ) = exp[-S off ({SJ)], 



(3) 



II. 



MODEL AND METHOD 



where S c fi is the effective action calculated as 



In this section, we introduce the model and method. 
The model is given in Sec. Ill Al and the MC method is 
described in Sec. Ill Bl In Sec. Ill Cl we give the definitions 
of physical quantities that we used to elaborate the phase 
diagram and thermodynamic properties. 



A. Model 

We consider a single-band Kondo lattice model on a 
triangular lattice with localized Ising spin moments. The 
Hamiltonian is given by 



* = -*£( 



C i<j C i° 



H.c.) + J$>*S;. (I) 



S*({Si}) = -£log[l + exp{-/3(^({ ) 5,})-M)}](4) 

V 

Here, E„({Si}) are the energy eigenvalues for the config- 
uration {Si}, which are readily calculated by the exact 
diagonalization as it is a one-particle problem in a static 
potential. 

The calculations were conducted for the system sizes 
TV = 12 x 12, 15 x 15, 12 x 18, and 18 x 18 under the pe- 
riodic boundary conditions. Thermal averages of physi- 
cal quantities were calculated for typically 4300-9800 MC 
steps after 1700-5000 steps for thermalization. The re- 
sults are shown in the temperature range where the ac- 
ceptance ratio is roughly larger than 1%. We divide the 
MC measurements into five bins and estimate the statis- 
tical errors by the standard deviations among the bins. 



C. Physical quantities 

As we will see later, the model (Q]) exhibits phase 
transitions to various magnetic states including different 
types of three-sublattice orders: ferrimagnetic (FR) state 
[Fig.fjjb)] and partially disordered (PD) state [Fig.QJc)]. 
These magnetic states, in principle, are distinguishable 
by the spin structure factor for the Ising spins, 



5(q) = ^E^^) ex P( it i- r '. 



(5) 



where the braket denotes the thermal average in the 
grand canonical ensemble, and ry is the position vector 
from i to jth site. The PD order is signaled by peaks of 
5(q) at q = ±(2tt/3, -2tt/3), while the FR order devel- 
ops a peak at q = in addition to q = ±(27r/3, — 27r/3). 
No Bragg peaks develop in the KT state as it is a quasi- 
LRO. However, in finite-size calculations, it is difficult to 
distinguish these phases solely by the structure factor, as 
the correlation length in the KT state is divergent and 
easily exceeds the system size at low temperature. 

For distinguishing the FR, PD, and KT instabilities, it 
is helpful to use the pseudospin defined for each three-site 
unit cell: 



(6) 
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and its summation 



™ = ^£s> 



(7) 



where m is the index for the three-site unit cells, and 
(i,j,k) denote the three sites in the mth unit cell 
belonging to the sublattices (A,B,C), respectivel y 12 ' 13 . 
Then, the three-sublattice PD state [Fig. Qlc)] is char- 
acterized by a finite M = (M x ,M y ,M z ) parallel to 
(•y/372, l/v^2,0), (0, a/2, 0), or their threefold symmet- 
ric directions around the z-axis. On the other hand, the 
three-sublattice FR state [Fig. QJb)] IS characterized by 
a finite M along (1/2A V% 1/v^), (2^273,0, -1/V3), 
or their threefold symmetric directions around the z-axis. 
Hence, the two states are distinguished by the azimuth of 
M in the xy-plane as well as M z . In the MC calculations, 
we measure 



M xy = ((M 2 + M 2 )^ 2 ), 
M z = (\M Z \), 

and the corresponding susceptibilities, 



JY 



X*v = 7 f ;{(M z x +M v i )-M* y ) 



T 

N 

T 



((M z 2 ) - M 2 Z ). 



(8) 
(9) 



(10) 
(11) 



We also introduce the azimuth parameter of M defined 
by 



ip = M cos 60a/, 



(12) 



where <pM is the azimuth of M in the xy plane and 
M = |M 2 . The parameter ip has a negative value and 
ip — > — 1| for the perfect PD ordering, while it becomes 
positive and i/j — > 1 for the perfect FR ordering; ip = 
for both paramagnetic and KT phases in the thermody- 
namic limit N — > oo. 

In addition, we calculate the spin entropy to distin- 
guish the three-sublattice orderings. The spin entropy 
per site is defined by 



5 ( T ) = -^E p ({^» lo g p ({^)' 



(13) 



{s,} 



where P({Si}) is the probability for spin configuration 
{Si} to be realized, given in Eq. ([3]). In the actual MC 
calculation, instead of directly calculating Eq. (TIB")) , S is 
evaluated by calculating its temperature derivative 



as(T) 



i 



dT NT 2 
and integrating it as 



{(S oS H)-(S oB )(H)} 



(14) 



sen = r 

Jo 



dS(T) 
dT 



dT = log2- 



dS{T) 
dT 



dT. (15) 



In Eq. (TT4")) . S e g is the effective action in Eq. (U). In the 
following calculations, we set the cutoff T = 1 for the 
upper limit of the last integral in Eq. (|15j) . 

On the other hand, in order to identify the two- 
sublattice stripe order [Fig.QJa)], we calculate the order 
parameter 



M s1 



E 



N 



1/2 



(16) 



and its susceptibility Xstr- Here, the sum is taken for the 
characteristic wave vectors of the stripe orders running in 
three different directions, q* tI . = (ir, 0) and (±571", ^7r). 
We also examine the thermodynamic behavior of elec- 
tronic states for itinerant electrons. There, we computed 
the charge modulation defined by 



"CO = 



|JV(q£o)| 



1/2 



(17) 



at q£.o = (— 27r/3, 27r/v / 3), which corresponds to the 
wave numbers for the three-sublattice orders. Here, iV(q) 
is the charge structure factor for itinerant electrons, 



N(q) 



— E(^ n j) ex p(iq- 



0. 



(18) 
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FIG. 2. (Color online). Mean-field band structure calcu- 
lated by Eq. (|19[) for the local magnetic field of PD type, 
(Aa, Ab, Ac) = (2, 0, —2). Each of the three bands shown is 
doubly degenerate, and there are totally six bands. The gray 
hexagon on the basal plane shows the first Brillouin zone for 
the magnetic supercell. 



III. 



MEAN-FIELD BAND STRUCTURE 



Before going to the MC results, we here discuss how 
one particle band structure is modulated by PD ordering 
in a mean-field picture. We consider a three-sublattice 
LRO state, in which the localized spins give a mean-field 
local magnetic field to itinerant electrons. Namely, we 
consider a mean-field Hamiltonian given by 
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(19) 



Here, three rows correspond to the different sublattices 
A, B, and C in the three-site unit cell; A a is a mean field 
given by J(S a ) (a = A,B,C). The sum is taken in the 
first Brillouin zone for the magnetic unit cell for three- 
sublattice order. Tk is the hopping term for itinerant 
electrons given by 



r k = -t[e lk 



.vg t 



. sdt 



(20) 



and a z a corresponds to the z component of itinerant elec- 
tron spin in each sublattice a. 

The band structure for a FR order, (Aa, Ab, Ac) = 
(A, A,— A), was recently studied by the authors^. 
There, it was reported that the electronic structure in 
the FR order is semimetallic with forming Dirac nodes at 
the electron filling n = j^ Z)i<t( c L c ^) = 1/3 10r J > *• 

Here, we discuss the band structure for the PD case, 
(A A ,A B ,A C ) = (A,0,-A). The band structure for 
A = 2 is shown in Fig. [21 In this case, all three bands 
shown in the figure arc doubly degenerate and there are 
six bands in total. The first Brillouin zone is shown by 
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FIG. 3. (Color online). Mean-field band structure along 
the symmetric lines in the local magnetic field of PD type, 

(Aa,Ab,A c ) = (A,0,-A): (a) A = 1/3, (b) A = 2/3, and 
(c) A = 2. The dashed horizontal lines indicate the Fermi 
level for n — 1/3. 



the gray shade in the bottom surface. The result shows 
the presence of an energy gap at the Fermi level cor- 
responding to n = 1/3, that opens between the lowest 
energy band and the middle band [see also Fig. [3fc)]. 

We next look into the conditions for the energy gap 
formation in the mean-field PD band. Figure [3] shows 
the results of band structure while varying A. The results 
are plotted along the symmetric line in the Brillouin zone 
shown in the bottom surface in Fig. [2] For small A, the 
system is metallic at n = 1/3, as shown in the case of 
A = 1/3 in Fig. E{a); both electron and hole pockets 
are present at the Fermi level. The pockets shrink as 
increasing A, and disappear at the same time at A = 2/3, 
as shown in Fig.^b). For larger A, an energy gap opens 
between the lowest and middle bands, corresponding to 
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FIG. 4. (Color online). A dependences of the mean-field 
energy gap and associated charge modulation nco at n — 1/3. 



n = 1/3, as stated above [Fig. Hfc)]. Hence, A c = 2/3 
is the critical point for the metal-insulator transition in 
this mean-field PD state. 

Figure [4] shows A depedences of the energy gap and 
associated charge modulation n co [Eq. ([T7|) ] at n = 1/3. 
The charge gap develops for A > 2/3 and monotonically 
increases, approaching asymptotically a A-linear form as 
A ^> t. The charge modulation is induced by the in- 
homogeneity of local potential; the local charge density 
at B sites (the site corresponds to paramagnetic sites) 
becomes dilute compared to those at A and C sites (the 
magnetically ordered sites). In the limit of A 3> t, nco 
approaches nco = l/\/l2 ~ 0.289. 

The results above suggest a stabilization mechanism of 
PD which is absent in the localized spin only model. In 
the previous studies on the Ising spin models^ ' 10 i 12 and 
an equivalent classical particle model^ on a triangular 
lattice, PD was shown to be unstable against thermal 
fluctuations and taken over by a KT state. In the case of 
our model, however, as the KT state lacks a long-range 
periodic magnetic structure, it is expected that the KT 
state does not open an energy gap in the electronic state 
of itinerant electrons. Therefore, in contrast to the case 
of localized spin only models, there is a chance for the 
current model to stabilize the PD state by the Slater 
mechanism, that is, by forming an energy gap at the 
Fermi level with folding the Brillouin zone under a peri- 
odic magnetic order. 

In addition, the formation of an energy gap for A > 
2/3 implies that, if the PD state is stabilized by the Slater 
mechanism, it should appear from a finite J, and not re- 
main stable down to J — > 0. This is in sharp contrast 
to magnetic ordering by the Ruderman-Kittel-Kasuya- 
Yosida (RKKY) interactional^; as the RKKY interac- 
tion is given by the second-order perturbation in terms 
of J/t, if the PD state is stabilized by the RKKY inter- 
action, it should appear for an infinitesimal J. Hence, 
the phase diagram in the small J region gives an idea on 
how the PD state is stabilized. We will discuss this point 
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FIG. 5. (Color online). Phase diagrams of the model (JlJ while 
varying n at (a) J = 1 and (b) J = 2. The symbols show 
phase boundaries for the four phases: stripe, partially disor- 
dered (PD), KT-like ("KT"), and ferrimagnetic (FR) phases. 
PS represents a phase separation. The lines are guides for the 
eyes. The strips at T — show the ground states obtained by 
comparing the energy of stripe and FR states. 



by showing the MC results while changing J in the next 
section. 



IV. 



MONTE CARLO SIMULATION 



In this section, we present the results of MC simu- 
lation introduced in Sec. Ill Bl We first show the finite- 
temperature phase diagrams in Sec. IIV A[ which include 
four magnetic phases: stripe, PD, FR, and KT-like 
states. The details of numerical data for the PD state 
are elaborated in Sec. IIVBI The results for stripe, KT- 
like, and FR states are discussed in Sec. IIV CI 



A. Phase diagrams 

Figure [5]Ja) shows the phase diagram around the elec- 
tron filling n = 1/3 at J = 1 obtained by MC calcula- 
tions. There are four dominant ordered phases — stripe, 
FR, PD, and KT-like phases, in addition to an elec- 
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FIG. 6. (Color online). Phase diagram of the model fTJ at 
n = 1/3 while varying J. The notations are common to those 
in Fig. [5] The boundary between PD and PS is difficult to 
determine by MC calculations, and supposed to be located at 
lower temperature than indicated by the gray arrows. 



tronic phase separation (PS). The strip at the bottom 
of the figure shows the ground state obtained by vari- 
ational calculation comparing the ground state energy 
of the stripe and FR states (the details of variational 
calculation is given in Appendix [A]) . For the relatively 
low filling of n < 0.29, the stripe order with period two 
[Fig. [TJa)] develops in the low temperature region. On 
the other hand, for the higher filling of n > 0.32, the sys- 
tem exhibits the three-sublatticc FR order at low tem- 
perature [Fig. [T|b)]. MC data for the stripe and FR 
orders will be discussed in Scc. lIVCl In addition to these 
two states, the numerical results show two intermediate- 
temperature states depending on the electron filling n. 
For 0.29 < n < 0.34, we identify the intermediate phase 
as the three-sublattice PD state [Fig. [T|c)] . The details 
will be discussed in Sec. IIVBI Meanwhile, for n > 0.34, 
we find KT-likc behavior similar to the one discussed in 
the Ising models^—, as presented in Sec. IIVC1 In these 
intermediate-temperature phases, the numerical data in- 
dicate a LRO for PD but a quasi-LRO in the KT-likc 
region. 

A similar phase diagram is obtained at J = 2, as shown 
in Fig. [5jb). In this case also, the PD phase emerges in 
the intermediate-temperature region. However, in con- 
trast to the case with J = 1 where PD is found widely 
above the FR state as well as PS, the PD phase domi- 
nantly appears above the PS region between the stripe 
and FR states. 

We also investigated the phase diagram of the model in 
Eq. (TTJ) while varying J . Figure [5] shows the numerically 
obtained phase diagram at n = 1/3. The result shows 
that the PD state is stable in a wide range of 0.8 < J < 
5.6. The transition temperature first rapidly increases as 
increasing J, while it turns to a gradual decrease after 
showing a peak at J ~ 2. 

An important observation in this constant-n phase di- 
agram is that the PD state docs not survive down to 



J — > 0, and it is taken over by the KT-like and FR phases 
in the small J region. The absence of PD state in the 
J — >• limit implies that the RKKY interaction in the 
second-order perturbation theory is insufficient in sta- 
bilizing the PD state. Moreover, the emergence of PD 
for J > J c =/= is consistently understood within the 
Slater mechanism discussed in Sec. Mlt the MC result of 
J c ~ 0.8 is in good accordance with the mean-field argu- 
ment of the critical value A c = 2/3. The result clearly 
indicates that a non-perturbative effect of itinerant elec- 
trons plays a crucial role in stabilizing the PD state. 

In the PD region in Fig. |6l our MC data do not show 
clear sign of further transition while decreasing temper- 
ature before the MC calculations become unstable. In 
the low temperature region, however, it becomes difficult 
to determine the chemical potential /i for n = 1/3. The 
lowest temperature of MC calculations are shown in the 
phase diagram by the gray downward arrows. On the 
other hand, the analysis of the ground state indicates 
that the ground state for J < 1.68 is the FR state, while 
the region for J > 1.68 is PS between the stripe and 
FR states. In addition, we observe the PS instability by 
carefully investigating the change of n as a function of 
(j, at J = 5.4 (see also Appendix |A")) . From these facts, 
we conclude that the PD for J > 1.68 is taken over by 
PS between the stripe and FR states. Since it is tedious 
to determine the PS boundary from fi-n plot for all the 
values of J, we merely plot the lowest temperature we 
reached in our constant-n calculations as the upper limit 
of temperature for the PS instability. 



B. 



Partial disorder 



Here, we present the details of MC data for identify- 
ing the PD state. Figure [7] shows T dependences of MC 
results for different J at n — 1/3. To fix n, we tuned 
[i for each temperature; the errors for n at each tem- 
perature are controlled within 0.001. Figure[7Jal) is the 
result for the pseudomoments M xy and M z at J = 1 [see 
the definitions in Eqs. flSJ) and ©, respectively]. M xy 
shows two anomalies while decreasing temperature at 



T 



PD) 



0.086(4) and T C (FR) = 0.019(2). The critical 
temperatures are determined by the peaks of the suscep- 
tibilities, Xxy, and Xz> as mentioned below. At X^ , 
M xy rapidly increases and approaches y/2 at lower tem- 
perature. In addition, it shows a kink at T c and fur- 
ther increase to 8/3 at lower temperature. Meanwhile, 
M z shows no anomaly at T c , while it shows a rapid 
increase to l/v3 at T c . Correspondingly, \xy and \z 
in Fig. [7|a2) also show divergent peaks increasing with 
the system size; peaks of Xxy appear at both Tc and 

Tc , while Xz shows a peak only at T c ■ These results 
signal the presence of two successive phase transitions 
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0.019(2). The error 
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FIG. 7. (Color online). MC results for (al)-(cl) M xy , M z , and tp, (a2)-(c2) \xy and Xz, and (a3)-(c3) S and its temperature 
derivative dS/dT at n — 1/3; (al)-(a3) J — 1, (bl)-(b3) J = 2, and (cl)-(c3) J = 4. The calculations were done for the system 
sizes N = 12 X 12, 12 x 18, and 18 x 18. <S is calculated from numerical integration of dS/dT by assuming S(T — 1) = log 2. 



the standard deviation of the MC data exceeds the dif- 
ference of expectation value from the peak value. The 
transition temperatures and error bars shown in Figs. [5] 
and [6] are given by this criterion. Meanwhile, most of the 
calculations in Fig. [5] were done by fixing fj, instead of n. 
Hence, we also give the error bars in terms of n, as n 
changes with T in a fixed y, calculation. 

To determine the nature of low temperature phases 
at n = 1/3, we also computed the azimuth parameter 
tp [Eq. (|12j) ] shown in Fig. Etal). While increasing the 
system sizes, ip apparently deviates from zero to a nega- 
tive value below T c , indicating that the intermediate 
phase for T C (FR) < T < T C (PD) has a PD type order. On 

( FRl 

the other hand, ip shows a sign change at T c , and 
rapidly increases to ip = 1 at lower temperature. This is 
a signature of the FR transition, which will be discussed 
in detail in Sec. ITVCl 

The emergence of PD is also seen in the results for the 
spin entropy S and its temperature derivative [Eqs. (fT5)) 
and (fl~4")) . respectively], as shown in Fig. [7Ja3). In the 

intermediate-temperature region for T c < T < T c , 
S appears to approach | log 2 as decreasing temperature, 
which is the value expected for the ideal PD state where 
one out of three spins in the magnetic unit cell remains 
paramagnetic. The remaining entropy is released rapidly 
at Tc and S — > at lower temperature due to the 
ordering of paramagnetic spins in the FR state. 

Similar phase transitions to the PD state are observed 
in the wide range of J, as shown in Figs.[7Jb) andjTJc) at 
J = 2 and J = 4, respectively. In these results, however, 
we could not confirm the presence of another phase tran- 
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FIG. 8. (Color online). MC results for S(q) along the q = 
(qx,0) line at T — 0.02. The calculations were done for the 
system size A*' = 18 x 18. 



sition at a lower temperature in the range of temperature 
we calculated, in contrast to the FR transition found in 
the case of J = 1 . As the PD state retains a finite S, it is 
unlikely that this phase survives to T — >• 0. Hence, it is 
presumably taken over by other ordered phases or PS at a 
lower temperature. As shown in Fig. |6l the ground state 
is deduced to be PS for the values of J in Figs. (7{b) and 
[7Jc). We, therefore, expect that the PD state is taken 
over by PS below T = 0.02 for J > 2. The situation 
is indicated by the gray arrows in the phase diagram in 
Fig. [?51 as discussed in Sec. IIVA1 

Another point to be noted is the systematic change in S 
in the PD state by changing J. While the result at J = 1 
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FIG. 9. (Color online). MC results for ip while varying n 
at T = 0.08 and J — 2. The calculations were done for the 
system sizes N = 12 x 12, 12 x 18, and 18 x 18. 



appears to show plateau like behavior at S ~ | log 2, the 
plateau value of 5 in the PD state decreases while in- 
creasing J, as shown in Figs. EJa3), [7Jb3), and[7Jc3). 
The decrease in S is presumably attributed to the de- 
velopment of spatial correlations between paramagnetic 
sites in the PD state; the ideal value S = | log 2 is for 
completely uncorrclated paramagnetic spins, and corre- 
lations between them reduces the entropy. Such devel- 
opment of correlatins are observed in the spin structure 
factor 5(q) defined in Eq. ([5]). Figure [5] shows a profile 
of S(q) calculated by MC simulation at T = 0.02. The 
peaks at q = (47r/3, 0) and (87r/3,0) indicates that the 
system is in a threc-sublattice ordered phase, while the 
absence of a sharp peak at q — (0, 0) indicates that there 
is no net magnetic moment; the result is consistent with 
PD order. When comparing the results at J = 2 and 
J = 4, the peak corresponding to the three-sublattice or- 
der gets sharper for J = 4, while the height of the peak 
of 5(q) is almost the same. This indicates that the PD 
order at J = 2 shows more spin fluctuations than that at 
J = 4, consistent with the trend of the plateau value of 
S. 

Thus far, we showed the results at n = 1/3. Next, 
we show how the PD evolves while changing n. Figure [9] 
shows the MC result of ip as a function of n at T = 0.08 
and J = 2. ifi becomes negative around n = 1/3 and 
takes the lowest value at n ~ 1/3. The data indicate that 
-0 is almost system size independent or rather slightly 
decreases as the system size increases in the finite range 
of n around n — 1/3. Hence, the PD state is stabilized 
not only at n = 1/3 but for a finite range of 0.31 < n < 
0.34 in the thermodynamic limit. The range well agrees 
with that for the PD phase estimated from the peak of 
susceptibilities shown in Fig.[S{b). 

With regard to the order of the PD transition, the PD 
transition in our MC results appears to be continuous, 
as shown in Fig. [7J However, it needs careful considera- 
tion, as we will discuss here. It is known that the Ising 
model on a triangular lattice with AF NN interactions 



is effectively described by a six-state model, in which 
the low-energy states with three up-up-down and three 
up-down-down configurations in the three-site unit cell 
are described by six-state variables. The PD state in our 
model also retains six low-energy states with different up- 
down-paramagnctic configurations, and hence, the tran- 
sition to PD is expected to be classified in the framework 
of six-state models. However, from the argument of du- 
ality properties, it is prohibited that the six-state mod- 
els exhibit a single second-order transition for changing 
temperature^. For instance, a two-dimensional six-state 
clock model shows two KT transitions at finite temper- 
ature, without exhibiting true LRO for T ^ 022i2i. On 
the other hand, a six-state Potts model shows a weak first 
order transition to LRO, in which the correlation length 
reaches the order of 1000 sites at the critical point2£ . In 
our PD case, the apparently second-order transition at 
T c is not expected to be a single one, but is always 
followed by another transition to FR or PS at a lower 
temperature. This appears not to violate the general ar- 
gument for the six-state models, although it is not clear 
to what extent the argument applies, as the electronic PS 
never takes place in the localized spin models. Hence, the 
PD transition can be of second order, as indicated in our 
numerical results. Of course, we cannot exclude the pos- 
sibility of a weak first order transition, similar to that 
of the Potts model. In this case, due to a long correla- 
tion length at the critical temperature, the system sizes 
used in our calculations arc likely to be insufficient to 
distinguish the first order transition from second order 
one. 



C. Other magnetic orders 

Figure [10] presents the results for the relatively low 
filling where the stripe order is stabilized at low temper- 
ature. Figure UOf a) shows the order parameter for the 
stripe order, M str [Eq. (|16p], and Fig. fTUTb 1 ) shows the 
corresponding susceptibility x str at J = 2 and n = 0.27. 
A phase transition to the stripe phase is signaled by a 
rapid increase of M str and corresponding peak of Xstr! 
we determine the transition temperature T c by the 
peak temperature of Xstr for each system size, and plot 
them in the phase diagram in Fig. [5Ja). The error bars 



arc estimated in a similar manner to the case of T, 



(PD) 



and T, 



(FR) 



We also show the system-size extrapolation 



of TJ str) in the inset of Fig. fTDTb). Although the data 
are rather scattered, we fit them by f(N) = a + b/N c 
with fitting parameters a, 6, and c. The extrapolation 
clearly shows that the phase transition takes place at a 
finite temperature, as expected for the two-dimensional 
Ising order. 

The stripe ordered phase is a peculiar magnetic state, 
in which the sixfold rotational symmetry of the lattice is 
spontaneously broken and reduced to twofold. Due to the 
symmetry breaking, the transport property is expected to 
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FIG. 10. (Color online). MC results for (a) M atr and (b) its 
susceptibility Xstr at J = 2 and n — 0.27. The inset in (b) 
shows Tc for different sizes and the solid line is the extrap- 
olation which gives T c = 0.051(13). The calculations were 
done for the system sizes N = 12 x 12, 14 x 14, 12 x 18, 16 x 16, 
and 18 x 18. 



show strong spatial anisotropy; e.g., the longitudinal con- 
ductivity will be large in the direction along the stripes, 
while suppressed in the perpendicular direction. This is 
an interesting topic on the control of transport by mag- 
netism and vice versa. 

Figure [TT1 shows the results for the relatively high fill- 
ing where the low temperature phase is FR, at n — 0.38 
and J = 2. The data indicate two successive tran- 
sitions signaled by the peaks in Xxy and Xz at differ- 
ent temperature. The peak of Xz corresponding to the 
increase of M z signals the phase transition to the FR 



phase at T< 



(FR) 



0.098(4). At the same time, ip be- 

,(FR) 



comes finite below T c , and approaches 1, as expected 
for the FR ordering. Similar behavior was observed at 



T C (FR) = 0.019(2) in Figs. [Hal) and[3;a2). On the other 
hand, at a higher Tkt = 0.146(4), only M xy changes 
rapidly, and correspondingly, Xxy shows a peak. M xy , 
however, shows a noticeable system-size dependence even 
below Tkt, in contrast with the results below T c 
Similar behavior was observed in the KT transition in 
Ising spin system a 12 ' 13 . 

On the other hand, ip does not show an anomaly at 

(FFM 

Tkt , while it shows a sharp rise around T c , as shown 
in Fig. fTTTa). The value of -0 extrapolated to large N 
converges to zero in the intermediate-temperature range. 
Figure [T2l shows the extrapolation of ip for TV — >• oo. The 
results indicate that, ip remains to be zero at N — > oo 
for T > 0.104, which is far below T KT = 0.146(4). On 
the other hand, the extrapolated value becomes finite 
for T < 0.104, reflecting the FR order; the transition 
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FIG. 11. (Color online). MC results for (a) M xy , M z , and 
i*: 3 ) Xxy and Xz, and (c) S and its temperature derivative 
dS/dT at 7i = 0.38 and J = 2. The calculations were done 
for the system sizes N - 12 x 12, 12 x 18, and 18 x 18. 




temperature is estimated as T, 



(FR) 



0.102(2), which is 
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FIG. 12. (Color online). Extrapolation of tp to N — > oo at 
different temperatures. The solid lines for T < 0.104 is the 
linear fitting of data. 



in accordance with T C (FR) = 0.098(4). 

The results above indicate that there is no sixfold 
symmetry breaking in M xy at Tkt, as seen in the KT 
phase in the Ising spin models^ 2 -. Hence, we consider 
that the higher-temperature transition at Tkt is of KT 
type. Namely, the system exhibits two successive transi- 
tions from the paramagnetic phase to the KT-like phase 
at Tkt, and the KT-like phase to the low-temperature 

(FR^ 

FR phase at T c • Here, we call the intermediate- 
temperature phase the KT-like phase, as it is difficult 
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FIG. 13. (Color online). MC results for the real-space spin 
correlation function C(r) at J = 2 and n = 0.38. The results 
are shown only for the sites with C(r) > 0. The calculations 
were done for the system size N — 18 X 18. 



to confirm cither the KT universality class by critical be- 
havior or the quasi-LRO behavior within the system sizes 
we reached, as seen below. 

The signature of two successive transitions is also ob- 
served in the real-space spin correlation function C(r). 
Here C{r) is the averaged correlations between the Ising 
spins in distance r, defined by 



c w=Ejv^ ( ^ (|1 



r), 



(21) 



where N p (r) = Xii7^(l r yl — r ) IS the number of spin 
pairs with distance r, and S(x) is the delta function. 
The MC data while varying temperature are shown in 
Fig. [13] Although the results are not conclusive due to 
the limitation on accessible system sizes, they appear to 
be consistent with the two transitions discussed above. 
For T < T C (FR) = 0.098(4), the spin correlation appears 
to approach constant for large distance, well correspond- 
ing to the FR LRO developed in this low temperature 
region. On the other hand, for T > T KT = 0.146(4), it 
becomes concave downward with a steep decrease with 
respect to the distance, which reflects an exponential de- 
cay in the high temperature paramagnetic state. In the 
intermediate region for T c < T < Tkt, the spin cor- 
relation also decays with increasing distance. The decay, 
however, is much slower and appears to obey an asymp- 
totic power law, which is characteristic to the quasi-LRO 
in the KT state. In principle, the critical exponents can 
be estimated from the asymptotic power-law behavior, 
but it is difficult to be conclusive in the current system 
sizes. 



V. ELECTRONIC STRUCTURE OF 
PARTIALLY DISORDERED STATE 

In the previous section, we discussed the thermody- 
namic behavior of the localized spin degree of freedom, 
with emphasis on the emergence of peculiar PD state. 
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FIG. 14. (Color online). MC results for nco at q = 
(2tt/3, -27r/3) atn= 1/3 and (a) J = 1, (b) J = 2, and 
(c) J = 4. The calculations were done for the system sizes 
N = 12 x 12, 12 x 18, and 18 x 18. 



In this section, we focus on the behavior in the charge 
degree of freedom of itinerant electrons in the PD phase. 

Figure [Ml shows temperature dependence of the charge 
modulation nco [Eq. (fTT)) ] at n = 1/3 for different J. 
Figure Q3Ia) is the result at J = 1 for different sys- 
tem sizes. The result shows an increase of nco below 
T ~ ri PD) = 0.086(4), indicating that the PD state is 
accompanied by charge modulation with period three. 
Similar onsets of charge modulation at T c are ob- 
served for larger J, as shown in Figs. Qjjb) and [T4Tc) ; 
the amplitude of the modulation in the PD phase in- 
creases monotonically as J increases. The magnitude of 
the charge modulation is in the same order compared 
to the mean-field result in Fig. [4j, while the growth is 
considerably suppressed by a factor of two to four. 

We next look into the electronic density of states 
(DOS) at different temperature. Figure [T5] shows the 
results for DOS while varying temperature at J = 2 and 
n = 1/3. The Fermi level is set at e — 0. Here, DOS was 
calculated by counting the number of energy eigenvalues 
as the histogram with the energy interval of 0.0375. In 
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FIG. 15. (Color online). MC results for DOS of itinerant 
electrons at n = 1/3 and J = 2 for N = 18 x 18. The Fermi 
level is set at e = 0. The statistical errors are comparable to 
the width of the lines. 



the paramagnetic region for T > T C (PD) = 0.130(4), DOS 
is featureless near the Fermi level. On the other hand, 
below Tc , an energy gap develops at the Fermi level 
for n = 1/3. The result shows that the PD state is an in- 
sulator, which supports the scenario that PD is stabilized 
by the Slater mechanism described in Sec. Mil Similarly 
to the charge modulation, the energy gap in the MC re- 
sults is largely suppressed compared to that obtained by 
the mean-field analysis in Fig. 2] This appears to show 
the importance of appropriately taking into account of 
thermal fluctuations. 



VI. 



SUMMARY 



To summarize, by a combined analysis of the mean- 
field type calculation and Monte Carlo simulation, we 
have investigated the origin of the partial disorder in the 
Ising-spin Kondo lattice model in a two-dimensional tri- 
angular lattice. In the mean-field type calculation, we 
have clarified that a local magnetic field of the partial 
disorder type induces a metal-insulator transition at 1/3 
filling at a critical value of the field. The result suggests 
that the three-sublattice partial disorder can give rise to 
an energy gap, and therefore, it has a chance to be stabi- 
lized through the Slater mechanism. On the other hand, 
in the Monte Carlo simulation, we have provided convinc- 
ing numerical results on the emergence of partial disor- 
der at finite temperatures where the stripe phase and 
the ferrimagnetic order compete with each other. The 
Monte Carlo result shows that the partially disordered 
state appears above a nonzero value of the spin-charge 
coupling, and that it is insulating and accompanied by 
charge disproportionation. The nonzero critical value of 
the spin-charge coupling and the opening of the charge 
gap arc both qualitatively consistent with the mean-field 
analysis. The results indicate that the partial disorder is 
stabilized by the Slater mechanism which is characteristic 




FIG. 16. (Color online), (a) The grand potetial Q and (b) 
electron filling n with respect to the chemical potential \x, 
numerically calculated by exactly diagonalizing the one-body 
Hamiltonian for itinerant electrons. The results are obtained 
at J = 2 with A s = 24 x 24 site superlattice of N = 12 x 12 
site unit cells. The strip at the left side of (b) shows the 
ground state at the corresponding filling. 



to itinerant magnets. Our results not only clarify the new 
mechanism of partial disorder in two dimensions but also 
pave the way for understanding of the interesting physics 
related to the peculiar coexistence of magnetic order and 
paramagnetic moments in itinerant electron systems. 

An interesting extension of the current work would be 
to consider the effect of quantum fluctuation of localized 
spins. In our result, the partial disorder remains stable 
down to very low temperature, implying that the para- 
magnetic spins are largely fluctuating and sensitive to 
perturbations at low temperatures. Hence, an interesting 
possibility is that, by including quantum fluctuations, the 
partial disorder is further stabilized and remains stable 
even in the ground state. Indeed, a similar partial disor- 
der was found in the ground state of the Kondo lattice 
model with quantum spins at half fillings*. Therefore, it 
is intriguing to examine the effect of quantum fluctua- 
tions on the present model with Ising spins. However, it 
is not straightforwardly calculated by the present Monte 
Carlo method. The interesting problem is left for future 
study. 
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Appendix A: Phase separation 

In this appendix, we present how to identify the PS 
region. First, we show the method we used to deter- 
mine the ground state phase diagram shown in Figs. [5] 
and [5] The ground state is obtained by variational cal- 
culations, i.e., by comparing the grand potential per site, 
il = (H)/N — [in, where \x is the chemical potential and 



n is the electron filling. Here, we compare $7 calculated 
for the magnetically ordered states, stripe and FR, which 
appear in the MC simulation at low temperature in the 
present parameter regions. The procedure is shown in 
Fig.ll6lat J = 2. FiguresfTfJta) andfTfJlb) show the results 
of 17 and n, respectively, calculated for stripe and FR or- 
ders. For n < —1.87 (yU > —1.87), f2 for the stripe order is 
lower (higher) than that for the FR order, indicating that 
the stripe (FR) state is the ground state in this region. 
At the critical value of /i ~ —1.87, the electron filling 
for the two states take different values, n ~ 0.301 in the 
stripe state and n ~ 0.334 in the FR state, as shown in 
Fig. fTBlf b) . This indicates that n changes discontinuously 
from n ~ 0.301 to n ~ 0.334 at the transition between 
the stripe and FR states. In other words, the system is 
unstable in the region of 0.301 < n < 0.334 against PS 
between the two states; the range of n is identified as the 
electronic PS. The PS regions in Fig. [5] are determined 
in this manner. Meanwhile, the PS region at n = 1/3 in 
Fig. [5] is identified by the similar calculations by changing 
J. 

Next, we describe how the PS region is determined at 
finite temperature in the MC calculation. In the MC sim- 
ulation using the grand canonical ensemble, PS is char- 
acterized by a sudden jump of n while sweeping [i. Fig- 
ure 1 171 shows a typical MC result for n as a function of //. 
The result at T = 0.048 shows a smooth change of n in 
the entire region of fi in the figure. On the other hand, 
the results at T = 0.040 and 0.044 show a sudden change 

from n ~ 0.290 to 0.315 at fj, 1.996. We roughly 

estimate the PS region by the values of n at the both 
ends of the jump. The results arc plotted in the phase 
diagrams in Fig. [5] The range of PS slightly depends on 
the system size, and hence, we plot the threshold values 
of n for each system size in the phase diagram. 
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